!! **********************************************************************
!!  Set of subroutines to calculate:
!!   -The velocity field of a thin sheet (subroutine VELOCITYFIELD)
!!      from the boundary conditions (subroutine READ_BC)
!!   -Thickening of the sheet (subroutine THICKEN)
!! **********************************************************************

      SUBROUTINE VELOCITYFIELD (TEMPSMa, Dt, AX, BY, n, m, nn,
     +			vissup, visinf, visTer, visco, nincogn,
     +			tallmax, alfa, nitermax, average_Szz, u, v,
     +			BCfile, nbanda, ifixv)

CC  TEMPSMa[Ma]	If != 0 and nitermax>1 then uses the previous velocity field.
CC  Dt [s]   	Time increment, if the amount of deformation is too 
CC			large then Dt will be reduced to fit 'quotadef'
CC  AX [m]	Domain length in x
CC  BY [m]	Domain length in y
CC  n, m	Number of nodes (in x,y directions) minus 1: column index runs from 0 (west) to n (east)
CC  nn		Number of nodes (n+1)*(m+1)
CC  nincogn	2*nn
CC  nbanda	4(n+1)+7 
!!  nitermax	max. number of iterations for convergence between viscosity and velocity (if exceeded => divergence).
!!			If 0 =>  Lineal problem: viscosity doesn't depend on strain rate.
!!			If 1 =>  No lineal problem (but linear equation on each time step): 
!!				 Viscosity depends on previous velocity. No iteration to find a consistent velocity and viscosity.
!!			If >1 => No lineal problem: Iterating method recalculating the viscosity with the new velocity.
CC  tallmax	Criterium of convergence. If SUM(Dvel/vel) < tallmax => convergence. (p.e., =0.01).
CC  alfa	when nitermax>1 => v[t+Dt]=alfa*v[t+Dt]+(1-alfa)*v[t]    0 explicit;  1 implicit;  0.5 recommended
CC  visTer [Pa]	Array of thermal term of viscosity: 
CC			1/2 * stregth[Pa*m] / layer.thickness[m]
CC  visco [Pa*s] Array of effective viscosity of the layer: visTer / ref.strain.rate
CC			Returns the calculated viscosity (a function of the new strain rate).
CC  vissup [Pa*s]  Imposed upper limit for calculated viscosity 'visco'. (~10^25 Pa*s)
CC  visinf [Pa*s]  Imposed lower limit for calculated viscosity 'visco'. (~10^22 Pa*s)
CC  average_Szz [Pa=N/m2] Integral of (P(z)*dz)/thickness_layer between the surface (z=topo) and Zcomp.
CC                      Where P(z) is the vertical stress (integral of rho(z)*g*dz), 
CC			Zcomp is the depth of compensation, and thickness_layer=Zcomp+elevation. 
CC			The independent term of the thin_sheet equilibrium equation is the lateral gradient of -average_Szz.
!!			Example: in a thickenned crust and in a ridge (thinner lithosphere) this average_Szz is lower.
!!			The velocity field goes from the lowest to the highest average_Szz.
!!  u, v [m/s]	Arrays of velocity in x, y directions in the previous time step.
!!			Sorted from west to east and south to north.
!!			Returns the new calculated velocity.
!!  BCfile [char*84]  Boundary conditions file name 
! ifixv=1 => keep in a file the points inside body 1 => velocity fixed in these points
!		also in subrutine VISCOSITAT

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!      INTEGER nincogn,nbanda,Li
      CHARACTER*84 TITOL_BC,BCfile
      PARAMETER (pi=3.1415926535897932D0, FACTEMP=3.1536D7,
     +		FACVEL=3.1536D10, quotadef=8.D-2)
!     +		nincogn_SD=40482, nbanda_SD=475)			!! nincogn_SD=2*(n+1)*(m+1), nbanda_SD=4(n+1)+7
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: visold,velold,b
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: a
      DIMENSION u(nn), v(nn),average_Szz(nn),
     +		visTer(nn), visco(nn), vel(nincogn)
!     +		a(nincogn_SD,nbanda_SD), b(nincogn)

      ALLOCATE(velold(nincogn),STAT=istat1) 
      ALLOCATE(visold(nn),STAT=istat2) 
!	velold=0.D0
!	visold=0.D0
      IF (istat1>0.OR.istat2>0) 
     +      STOP " Error allocation SUBR. VELOCITYFIELD" 
      ALLOCATE(b(nincogn),STAT=istat3) 
      IF (istat1>0.OR.istat2>0.OR.istat3>0) 
     +      STOP " Error allocation SUBR. VELOCITYFIELD" 
      
      ALLOCATE(a(nincogn,nbanda),STAT=istat) 
      IF (istat>0) STOP " Error allocation a(:,:)" 
        
!       IF (nincogn_SD/=nincogn .OR. nbanda_SD/=nbanda)
!     +      WRITE(6,"(3X,'Matrix a not well dimensioned thin_sheet.f',
!     +      ' nincogn_SD/=',I6,' or nbanda_SD/=',I5/)") nincogn,nbanda


      tall=0.D0
      Dx=AX/n
      Dy=BY/m
      
      visold=visco
      vel=0
      kvel=1
      DO kixy=1,nn
          velold(kvel)=u(kixy)
          velold(kvel+1)=v(kixy)
          kvel=kvel+2
      END DO
CC *************** CONSTRUCCIO DE LA MATRIU *********************
      Li=2*(n+1)+3
      Ls=Li
      Ld=Li+1
      kp=Ld+2*(n+1) 
      kn=Ld-2*(n+1)
       
      niter=0
 111   CONTINUE 
       niter=niter+1           
       a=0.D0
       b=0.D0
CC  ----------  Condicions de Contorn --------------------
      CALL READ_BC (TEMPSMa,niter,n,m,nn,nincogn,nbanda,Dx,Dy,visco,
     +			a,b,TITOL_BC, BCfile)

C  ------------------------------------------------------
C   PUNTS INTERIORS
       DP_max=-1.D100
       DP_min=1.D100
       Dvis_max=-1.D100
       Dvis_min=1.D100
      Point_iy: DO iy=1,m-1
         Point_ix: DO ix=1,n-1
               k=ix+1+iy*(n+1)
                GLv=visco(k)
                G1=visco(k+1)
                G2=visco(k-1)
                G3=visco(k+n+1)
                G4=visco(k-n-1)
C  ---------------------------------------------------------------------
               GLvx=(G1-G2)/(2.D0*Dx)
               GLvy=(G3-G4)/(2.D0*Dy)
	       Dvis_max=MAX(Dvis_max,GLvx,GLvy)
	       Dvis_min=MIN(Dvis_min,GLvx,GLvy)
	       DP_x=(average_Szz(k+1)-average_Szz(k-1))/
     +	       		(2.D0*Dx)
	       DP_y=(average_Szz(k+n+1)-average_Szz(k-n-1))/
     +	       		(2.D0*Dy)
               DP_max=MAX(DP_max,DP_x,DP_y)
               DP_min=MIN(DP_min,DP_x,DP_y)
	       Tx=-DP_x
	       Ty=-DP_y
               l=2*ix+1+2*iy*(n+1)
               DIAGONAL1=-(8.D0*GLv)/(Dx*Dx)-(2.D0*GLv)/(Dy*Dy)
               a(l,Ld+2)=((4.D0*GLv)/(Dx*Dx)+(2.D0*GLvx)/Dx)/DIAGONAL1
               a(l,kp)=(GLv/(Dy*Dy)+GLvy/(2.D0*Dy))/DIAGONAL1
               a(l,Ld)=1.D0
               a(l,kn)=(GLv/(Dy*Dy)-GLvy/(2.D0*Dy))/DIAGONAL1
               a(l,Ld-2)=((4.D0*GLv)/(Dx*Dx)-(2.D0*GLvx)/Dx)/DIAGONAL1
               a(l,kp+3)=((3.D0*GLv)/(4.D0*Dx*Dy))/DIAGONAL1
               a(l,kn+3)=-a(l,kp+3)
               a(l,kp-1)=-a(l,kp+3)
               a(l,kn-1)=a(l,kp+3)
               a(l,kp+1)=(GLvx/Dy)/DIAGONAL1
               a(l,kn+1)=-a(l,kp+1)
               a(l,Ld+3)=(GLvy/(2.D0*Dx))/DIAGONAL1
               a(l,Ld-1)=-a(l,Ld+3)
               b(l)=Tx/DIAGONAL1
               l=2*ix+2+2*iy*(n+1)
               DIAGONAL2=-2.D0*GLv*(1.D0/(Dx*Dx)+4.D0/(Dy*Dy))
               a(l,Ld+2)=(GLv/(Dx*Dx)+GLvx/(2.D0*Dx))/DIAGONAL2
               a(l,kp)=((4.D0*GLv)/(Dy*Dy)+(2.D0*GLvy)/Dy)/DIAGONAL2
               a(l,Ld)=1.D0
               a(l,kn)=((4.D0*GLv)/(Dy*Dy)-(2.D0*GLvy)/Dy)/DIAGONAL2
               a(l,Ld-2)=(GLv/(Dx*Dx)-GLvx/(2.D0*Dx))/DIAGONAL2
               a(l,kp+1)=((3.D0*GLv)/(4.D0*Dx*Dy))/DIAGONAL2
               a(l,kn+1)=-a(l,kp+1)
               a(l,kp-3)=-a(l,kp+1)
               a(l,kn-3)=a(l,kp+1)
               a(l,Ld+1)=(GLvy/Dx)/DIAGONAL2
               a(l,Ld-3)=-a(l,Ld+1)
               a(l,kp-1)=(GLvx/(2.D0*Dy))/DIAGONAL2
               a(l,kn-1)=-a(l,kp-1)
               b(l)=Ty/DIAGONAL2
         END DO Point_ix 
      END DO Point_iy 

!! Fixing the velocity to the nodes of the file: 'P_fixVelocity.res'
!      IF(ifixv==1 .AND. TEMPSMa>=7.0 .AND. TEMPSMa<=38.0) THEN		!! Rif cap a l'Oest
!	   vx_fix=-8.D0		!! x velocity  [mm/yr]
!	   vy_fix=0.D0		!! y velocity [mm/yr]
      IF(ifixv==1) THEN		!! India indenter
	   vx_fix=0.D0		!! x velocity  [mm/yr]
	   vy_fix=50.D0		!! y velocity [mm/yr]
	   OPEN(17,FILE='P_fixVelocity.res',STATUS='OLD',ACTION='READ')
	       npoints=0
	       DO ip=1,nn
		   READ(17,*,END=440) ix,iy
		   npoints=npoints+1
		   leq1=2*ix+1+2*iy*(n+1)
		   leq2=leq1+1
		   DO icol=1,nbanda
		      a(leq1,icol)=0.D0
		      a(leq2,icol)=0.D0
		   END DO
		   a(leq1,Ld)=1.D0
		   a(leq2,Ld)=1.D0
		   b(leq1)=vx_fix/FACVEL
		   b(leq2)=vy_fix/FACVEL
	       END DO
 440	   WRITE(6,"(3X,'Number of interior points with the velocity ,'
     +			'fixed:',I5)") npoints
	   CLOSE(17)
      END IF

      CALL sistbanda(a,b,nincogn,nbanda,Li,vel)

!      IF((TEMPSMa==0.D0.AND.niter/=1).OR.(TEMPSMa/=0.D0.AND.nitermax>1))
      IF(TEMPSMa/=0.D0)
     +		vel=alfa*vel+(1.D0-alfa)*velold

      kvel=1
      DO iy=0,m
	 DO ix=0,n
	     kxy=ix+1+iy*(n+1)
	     u(kxy)=vel(kvel)
	     v(kxy)=vel(kvel+1)
	     kvel=kvel+2
         END DO
      END DO

       epmig=0.D0
       cepmig=0.D0
       epmax=0.D0
       NVISSUP=0
       NVISINF=0
      IF(nitermax==0) THEN
            PRINT*,'   No Iteration. The viscosity is the same'
              DO iy=1,m-1
                  DO ix=1,n-1
		      epeffec=effective_strainrate (ix, iy, m, n, nn,
     +							Dx, Dy, u, v)
		      epmax=MAX(ABS(epeffec),epmax)
	    	      epmig=epmig+ABS(epeffec)
		      cepmig=cepmig+1
                  END DO  
              END DO  
	    epmig=epmig/cepmig
            GOTO 999
      ENDIF
      IF(niter==1) WRITE(6,"(/4X,'Iteration',31X,
     +		'__ nodes limited for viscosity __'/
     +		4X'number',5X,'mean Dv',6X,'mean Dvis',7X,
     +		1P,G9.2,' Pa.s',4X,1P,G9.2,' Pa.s')") visinf,vissup

C ****************  TROBO LA NOVA VISCOSITAT *********************
      Vis_iy: DO iy=1,m-1
         Vis_ix: DO ix=1,n-1
		kxy=ix+1+iy*(n+1)
		epeffec=effective_strainrate (ix, iy, m, n, nn,
     +							Dx, Dy, u, v)
		epmig=epmig+ABS(epeffec)
		cepmig=cepmig+1
		epmax=MAX(ABS(epeffec),epmax)
		visnova=ABS(visTer(kxy)/epeffec)
		visco(kxy)=alfa*visnova+(1.D0-alfa)*visold(kxy)
		IF(visco(kxy).GT.vissup) THEN
			visco(kxy)=vissup
			NVISSUP=NVISSUP+1
		ENDIF
		IF(visco(kxy).LT.visinf) THEN
			visco(kxy)=visinf
			NVISINF=NVISINF+1
		ENDIF
         END DO Vis_ix
      END DO Vis_iy
      !!    FIXO LA VISCOSITAT A LES VORES 
           DO ix=1,n-1
                 kxys=ix+1
                 kxyn=ix+1+m*(n+1)
                 visco(kxys)=visco(kxys+n+1)
                 visco(kxyn)=visco(kxyn-n-1)
           END DO
           DO iy=0,m
                 kxyw=1+iy*(n+1)
                 kxye=n+1+iy*(n+1)
                 visco(kxyw)=visco(kxyw+1)
                 visco(kxye)=visco(kxye-1)
           END DO

      epmig=epmig/cepmig
      Dumig=0.D0
!      Dumax=0.D0
      velmax=0.D0
      Dvismig=0.D0
      Dvismax=0.D0
      NPU=0
      DO iy=1,m-1
         DO ix=1,n-1
              kd=ix+1+iy*(n+1)
              k=2*ix+1+2*iy*(n+1)
              velmod=DSQRT(vel(k)*vel(k)+vel(k+1)*vel(k+1))
              veloldm=DSQRT(velold(k)*velold(k)+velold(k+1)*velold(k+1))
              Du=(ABS(velmod-veloldm))/velmod  
              Dumig=Dumig+Du  
!              Dumax=MAX(Du,Dumax)
              Dvis=(ABS(visold(kd)-visco(kd)))/visco(kd)
              Dvismig=Dvismig+Dvis
              Dvismax=MAX(Dvis,Dvismax) 
              velmax=MAX(velmod,velmax) 
	      NPU=NPU+1
         END DO
      END DO
      Dvismig=Dvismig/NPU 
      Dumig=Dumig/NPU     
      vismax=VAL_MAX(visco,nn)
      IF(velmax.EQ.0.D0.OR.vismax.EQ.0.D0) THEN
              PRINT*,'             EL CAMP DE VELOCITATS ES NUL'
              PRINT*,' Terme maxim de la velocitat (velmax) =',velmax,
     +           '     Terme maxim de la viscositat (vismax) =',vismax
              GOTO 999 
      ENDIF
       tall=Dumig
       tallvis=Dvismig
       WRITE(6,"(4X,I3,5X,F10.5,4X,F10.5,12X,I4,14X,I4)") 
     +			niter,tall,tallvis,NVISINF,NVISSUP 
       if(niter.ge.nitermax) THEN    
             PRINT*,'   Maximum number of iterations'
             GOTO 999
       ENDIF

       IF(tall.GT.tallmax) then
	   visold=visco
           velold=vel
           GOTO 111
       ENDIF   
       WRITE(6,"(4X,'velocity converged')")     

999    CONTINUE 

      vismax=VAL_MAX(visco,nn)
      vismin=VAL_MIN(visco,nn)

      strainmx=epmax*Dt
      WRITE(6,61) vismin,vismax,epmig,epmax,epmig*Dt,strainmx
 61   FORMAT(
     +  4X,'Viscosity:            minim:',1P,G12.4,' Pa.s,   maxim:',
     +    	1P,G12.4,' Pa.s'/
     +	4X,'Effect. strain rate:  medium:',1P,G12.4,' s-1,   maximum:',
     +		1P,G12.4,' s-1'/
     +	4X,'Effect. strain:       medium:',1P,G12.4,'        maximum:',
     +		1P,G12.4)
     	 
          IF(strainmx.GT.quotadef) THEN	
               Dt=(quotadef/epmax)-5.D4
               Dtany=Dt/FACTEMP
               PRINT*,'TOO MUCH DEFORMATION -> DECREASE '
               PRINT*,'THE TIME INTERVAL,  Dt =',Dtany,' anys'
               strainmx=epmax*Dt
               PRINT*,'         NEW MAXIMUM STRAIN :',strainmx
          ENDIF

CC	Divideixo el vector vel (velocitat) en dos vectors:
CC	u (velocitat en x) i v (velocitat en y).
	kvel=1
      DO iy=0,m
	 DO ix=0,n
	     kxy=ix+1+iy*(n+1)
	     u(kxy)=vel(kvel)
	     v(kxy)=vel(kvel+1)
	     kvel=kvel+2
         END DO
      END DO
      
      IF (ALLOCATED(velold)) DEALLOCATE(velold,STAT=istat)       
      IF (ALLOCATED(visold)) DEALLOCATE(visold,STAT=istat)   
      IF (ALLOCATED(b)) DEALLOCATE(b,STAT=istat)
    
      RETURN
      END SUBROUTINE VELOCITYFIELD
       
CC ********************************************************************
CC ********************************************************************

      SUBROUTINE READ_BC (TEMPSMa,niter,n,m,nn,nincogn,nbanda,Dx,Dy,
     +				visco,a,b,TITOL_BC,BCfile)   

CC   Read the Boundary Conditions from the file: BC.in
CC   First line: Title
CC	ix,iy,ITBC,t1,t2 
CC	ITBC: Boundary Condition type.
CC	t1,t2: Boundary Condition 1 and 2.

CC  ITBC=1    velocity fixed [Vx(m/s), Vy(m/s)]			->  vel_x=t1  and  vel_y=t2
CC  ITBC=12   stress xx and xy fixed (Est, West free)		-> tau_xx=t1  and  tau_xy=t2
CC  ITBC=13   stress xx and yy fixed   		   		-> tau_xx=t1  and  tau_yy=t2
CC  ITBC=23   stress xy and yy fixed (North, South free)	-> tau_xy=t1  and  tau_yy=t2
CC  ITBC=4    free slip (vel normal=0, tau_xy=0)		-> don't use t1 and t2
CC  ITBC=5    velocity fixed [modul(mm/yr), azimut(degree)]	-> modul=t1  and  azimut=t2

CC  Stress free -> normal stress and xy zero.

    
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION a(nincogn,nbanda),b(nincogn),visco(nn)
      CHARACTER*84 TITOL_BC,BCfile
      PARAMETER (NBCmax=1000,FACVEL=3.1536D10,PI=3.1415926535897932D0)

      OPEN(4,file=BCfile)
      NBC=2*(n+1)+2*(m-1)
      Li=2*(n+1)+3
      Ls=Li
      Ld=Li+1
      kp=Ld+2*(n+1)
      kn=Ld-2*(n+1)
      READ(4,"(A80)") TITOL_BC
!     IF(TEMPSMa.EQ.0.AND.niter.EQ.1) WRITE (6,"(
!     	  WRITE (6,"(/2X,'READ BOUNDARY CONDITIONS: ',/4X,A80/
!     +    4X,'(ix, iy, BC type:   u,v / tau_1,tau_2 )')") TITOL_BC
      
       NBC_read=0
       DO 7 inBC=1,NBCmax
           READ(4,*,END=120) ix,iy,ITBC,t1,t2
!	    IF(TEMPSMa==0.AND.niter==1) 
!             	WRITE (6,"(1X,3I4,1P,2G14.6)") ix,iy,ITBC,t1,t2
           NBC_read=NBC_read+1
	   leq1=2*ix+1+2*iy*(n+1)
	   leq2=2*ix+2+2*iy*(n+1)
	   kxy=ix+1+iy*(n+1)
	   kp=Ld+2*(n+1)
           kn=Ld-2*(n+1)
CCC ====================================================
CCCC	ITBC=1 o 5  fixo (u,v)
	    IF(ITBC==1.OR.ITBC==5) THEN
	   	vel_x=t1
		vel_y=t2
		IF(ITBC==5) THEN	! read modul(mm/yr),azimut(degree)			
			vel_x=(t1/FACVEL)*DSIN(t2*PI/180.D0)
			vel_y=(t1/FACVEL)*DCOS(t2*PI/180.D0)
		ENDIF
             	a(leq1,Ld)=1.D0
             	b(leq1)=vel_x
             	a(leq2,Ld)=1.D0
             	b(leq2)=vel_y
		GOTO 99
	    ENDIF	   
CCC ====================================================
CCCC	ITBC=12  fixo:  tau_xx (eq.1) i  tau_xy (eq.2)
CCCC	ITBC=13  fixo:  tau_xx (eq.1) i  tau_yy (eq.2)
CCCC	ITBC=23  fixo:  tau_xy (eq.1) i  tau_yy (eq.2)
	    IF(ITBC.GT.11.AND.ITBC.LT.20) THEN
	   	tau_xx=t1
		GOTO 411
	    ENDIF
 412	    CONTINUE
	    IF(ITBC.EQ.12) THEN
		tau_xy=t2
		ieq=1
		GOTO 522
	    ENDIF
	    IF(ITBC.EQ.13) THEN
		tau_yy=t2
		GOTO 813
	    ENDIF
	    IF(ITBC.EQ.23) THEN
		tau_xy=t1
		ieq=0
		GOTO 522
	    ENDIF
 525	    CONTINUE
	    IF(ITBC.EQ.23) THEN
		tau_yy=t2
		GOTO 813
	    ENDIF
CCC ====================================================
CCCC	ITBC=4  free slip (v norm=0, dv(tang)/dx(tang) =0)
	    IF(ITBC.EQ.4) THEN
		IF(iy.EQ.0) THEN
             	    a(leq1,kp)=1.D0
             	    a(leq1,Ld)=-1.D0
             	    b(leq1)=0.D0
             	    a(leq2,Ld)=1.D0
		    b(leq2)=0.D0
		    GOTO 99
	    	ENDIF		
		IF(iy.EQ.m) THEN
             	    a(leq1,Ld)=1.D0
             	    a(leq1,kn)=-1.D0
             	    b(leq1)=0.D0
             	    a(leq2,Ld)=1.D0
		    b(leq2)=0.D0
		    GOTO 99
	    	ENDIF		
		IF(ix.EQ.0) THEN
             	    a(leq1,Ld)=1.D0
             	    b(leq1)=0.D0
             	    a(leq2,Ld+2)=1.D0
             	    a(leq2,Ld)=-1.D0
		    b(leq2)=0.D0
		    GOTO 99
	    	ENDIF		
		IF(ix.EQ.n) THEN
             	    a(leq1,Ld)=1.D0
             	    b(leq1)=0.D0
             	    a(leq2,Ld)=1.D0
             	    a(leq2,Ld-2)=-1.D0
		    b(leq2)=0.D0
		    GOTO 99
	    	ENDIF		
	    ENDIF		

CCCC---------------------------------------------------------
CCCC---------------------------------------------------------
CC  Condition   tau_xx:  Equation 1
 411		CONTINUE
		IF(ix.NE.0.AND.ix.NE.n) THEN	
		    a(leq1,Ld+2)=1.D0
		    a(leq1,Ld-2)=-1.D0
             	    b(leq1)=(Dx*tau_xx)/visco(kxy)
		ENDIF
		IF(ix.EQ.0) THEN	
		    a(leq1,Ld+2)=1.D0
		    a(leq1,Ld)=-1.D0
             	    b(leq1)=(Dx*tau_xx)/(2.D0*visco(kxy))
		ENDIF
		IF(ix.EQ.n) THEN	
		    a(leq1,Ld)=1.D0
		    a(leq1,Ld-2)=-1.D0
             	    b(leq1)=(Dx*tau_xx)/(2.D0*visco(kxy))
		ENDIF
		GOTO 412
CCCC---------------------------------------------------------
CCCC---------------------------------------------------------
CC  Condition	tau_xy:  Equation 1 (ieq=0),	  Equation 2 (ieq=1)
 522		CONTINUE
		leq=leq1+ieq
		IF(iy.NE.0.AND.iy.NE.m) THEN
		    a(leq,kp-ieq)=1.D0/(2.D0*Dy)
		    a(leq,kn-ieq)=-1.D0/(2.D0*Dy)
		ENDIF
		IF(iy.EQ.0) THEN
		    a(leq,kp-ieq)=1.D0/Dy
		    a(leq,Ld-ieq)=-1.D0/Dy
		ENDIF
		IF(iy.EQ.m) THEN
		    a(leq,Ld-ieq)=1.D0/Dy
		    a(leq,kn-ieq)=-1.D0/Dy
		ENDIF
		IF(ix.NE.0.AND.ix.NE.n) THEN
		    a(leq,Ld+3-ieq)=1.D0/(2.D0*Dx)
		    a(leq,Ld-1-ieq)=-1.D0/(2.D0*Dx)
		ENDIF		    
		IF(ix.EQ.0) THEN
		    a(leq,Ld+3-ieq)=1.D0/Dx
		    a(leq,Ld+1-ieq)=-1.D0/Dx
		ENDIF
		IF(ix.EQ.n) THEN
		    a(leq,Ld+1-ieq)=1.D0/Dx
		    a(leq,Ld-1-ieq)=-1.D0/Dx
		ENDIF
             	b(leq)=tau_xy/visco(kxy)
	    	IF(ieq.EQ.1) GOTO 99
	    	IF(ITBC.EQ.23) GOTO 525
CCCC---------------------------------------------------------
CCCC---------------------------------------------------------
CC  Condition   tau_yy:  Equation 2
 813		CONTINUE
		IF(iy.NE.0.AND.iy.NE.m) THEN	
		    a(leq2,kp)=1.D0
		    a(leq2,kn)=-1.D0
             	    b(leq2)=(Dy*tau_yy)/visco(kxy)
		ENDIF
		IF(iy.EQ.0) THEN	
		    a(leq2,kp)=1.D0
		    a(leq2,Ld)=-1.D0
             	    b(leq2)=(Dy*tau_yy)/(2.D0*visco(kxy))
		ENDIF
		IF(iy.EQ.m) THEN	
		    a(leq2,Ld)=1.D0
		    a(leq2,kn)=-1.D0
             	    b(leq2)=(Dy*tau_yy)/(2.D0*visco(kxy))
		ENDIF
	        GOTO 99
CCCC---------------------------------------------------------
		
 99	CONTINUE
 7      CONTINUE

 100  CONTINUE
      IF(NBC_read.EQ.0) PRINT*,' FILE NO READ. number of data:',NBC_read
      WRITE(6,"(/' Dataset contained more than ',I5,
     +           ' data; Check the dimensions.')") NBCmax
 120  CONTINUE
!      WRITE(6,"(/5X,'Reading of data completed:',I7,' data points'/ )")
!     +            NBC_read
 	
      IF(NBC_read.NE.NBC) THEN
           WRITE(6,"(' El numero de punts llegits:',I7,' no coincideix',
     +         'amb 2(n+1)+2(m-1)=',I7//'PROGRAMA ATURAT')")NBC_read,NBC
           STOP
      ENDIF
      CLOSE(4)

      RETURN
      END SUBROUTINE READ_BC

CC **********************************************************************
CC **********************************************************************

       SUBROUTINE THICKEN (Dt, n, m, nn, AX, BY, u, v, epuntzz, 
     +			       thickness, IBC_thicken)

!  New thickness of a layer with velocity field (u,v) after a time interval Dt(s)
!  (u,v) (m/s) : array of horizontal velocities.
!	d(thickness)/dt=
!		thickness*epuntzz-(u*d(thickness)/dx+v*d(thickness)/dy)

!      Amb els logaritmes dona problemes quan algun gruix es zero:
!	d[ln(thickness)]/dt=
!		epuntzz-(u*d[ln(thickness)]/dx+v*d[ln(thickness)]/dy)

! IBC_thicken=0 -> No temporal thickening variations on the boundaries,
!			thickness=thickness_old
! IBC_thicken!=0 -> No lateral variations of the thickening on the boundaries,
!			d(thickness)/dx=d(thickness)/dy=0


      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: thickness_old
      DIMENSION u(nn),v(nn), epuntzz(nn), thickness(nn)
      
      ALLOCATE(thickness_old(nn),STAT=istat) 
      IF (istat>0) STOP " Error allocation thickness_old(:)" 
      Dx=AX/n
      Dy=BY/m
      thickness_old(:)=thickness(:)

      Column: DO iy=1,m-1
         Row: DO ix=1,n-1
            kxy=ix+1+iy*(n+1)
            TERME1=thickness_old(kxy)*epuntzz(kxy)
            !sx=(thickness_old(kxy+1)-thickness_old(kxy-1))/(2.D0*Dx)		!! Derivada centrada
            !sy=(thickness_old(kxy+n+1)-thickness_old(kxy-n-1))/(2.D0*Dy)	!! Derivada centrada	    

	    sign_x=SIGN(1.D0,u(kxy))			!! Derivada segons el fluxe del fluid		
	    sign_y=SIGN(1.D0,v(kxy))			!! Derivada segons el fluxe del fluid
	    sx_sup=(thickness_old(kxy+1)-thickness_old(kxy))/Dx
            sx_inf=(thickness_old(kxy)-thickness_old(kxy-1))/Dx		
            sy_sup=(thickness_old(kxy+n+1)-thickness_old(kxy))/Dy
            sy_inf=(thickness_old(kxy)-thickness_old(kxy-n-1))/Dy	    	    
            sx=sx_inf*MAX(sign_x,0D0)-sx_sup*MIN(sign_x,0D0)
            sy=sy_inf*MAX(sign_y,0D0)-sy_sup*MIN(sign_y,0D0)
	    
!!           TERME2=0.D0 		! LAGRANGIAN sist. (grid moves with the material)
            TERME2=u(kxy)*sx+v(kxy)*sy	! EULERIAN sist. (grid fixed)
            thickness(kxy)=thickness_old(kxy)+Dt*(TERME1-TERME2)	    
	    IF(thickness(kxy).LT.0.0) thickness(kxy)=0.D0
         END DO Row
      END DO Column
!  Boundary Conditions     
      IF(IBC_thicken==0) THEN  ! ix=0,ix=n,iy=0,iy=m -> thickness=thickness_old
	  DO iy=1,m-1
             kxy0=1+iy*(n+1)
             kxyn=n+1+iy*(n+1)
             thickness(kxy0)=thickness_old(kxy0)
             thickness(kxyn)=thickness_old(kxyn)
          END DO
          DO ix=0,n
             kxy0=ix+1
             kxym=ix+1+m*(n+1)
             thickness(kxy0)=thickness_old(kxy0)
             thickness(kxym)=thickness_old(kxym)
          END DO
       ELSE   ! ix=0,ix=n,iy=0,iy=m -> d(thickness)/dx=d(thickness)/dy=0
          DO iy=1,m-1
             kxy0=1+iy*(n+1)
             kxyn=n+1+iy*(n+1)
             thickness(kxy0)=thickness(kxy0+1)
             thickness(kxyn)=thickness(kxyn-1)
          END DO
          DO ix=0,n
             kxy0=ix+1
             kxym=ix+1+m*(n+1)
             thickness(kxy0)=thickness(kxy0+n+1)
             thickness(kxym)=thickness(kxym-n-1)
          END DO
      ENDIF 
CC --------------------------------------------------------------------
!      CALL smooth_horizontal (m, n, nn, Dx, Dy, thickness)
      thickmin=VAL_MIN(thickness,nn)
      thickmax=VAL_MAX(thickness,nn)
      WRITE(6,"(4X,'Thickness:',12X,'minim:',F10.2,' m',9X,
     +        'maxim: ',F10.2,' m')") thickmin,thickmax
      
      IF (ALLOCATED(thickness_old)) DEALLOCATE(thickness_old,STAT=istat)     
      RETURN
      END SUBROUTINE THICKEN



